socp <-
function(A, b, c, kvec, type = rep('q',length(kvec)), delta0 = .75, sigma0 = .25, mu0 = 0.01, zero_tol = 1e-6, max_iter = 100, min_progress = zero_tol/100){
	#A is m x sum(kvec)
	#b is a sum(kvec) vector
	#kvec = c(k1, ..., kn)
	
	#Set return code
	code = -1
	#Codes are: 
	#0 - Full success, 
	#1 - Singularity occured,
	#2 - Lack of progress,
	#3 - Maximum number of iterations reached
	
	#Check the condition number of A
	condnum <- rcond(A, useInve = TRUE)
	if(condnum < zero_tol){
		cat("Reciprocal condition number of A is not within tolerance.\nrcond(A, useInv=TRUE) ==",condnum,"\n")
	}
	
	#Change linear constraints to SOC
	converted <- convert_l(kvec,type)
	kvec <- converted$kvec
	type <- converted$type
	
	#Dimensions
	n <- length(kvec)
	m <- dim(A)[1]
	k <- dim(A)[2]
	
	#Create the SOC index list
	inds <- make_inds(kvec, type)
	
	#Step 0
	delta <- delta0
	sigma <- sigma0
	mu <- rep(mu0,n)
	zbar <- c(rep(0,k + m), mu)
	x <- jordan_identity(kvec)
	y <- rep(0, m)
	s <- c - t(A)%*%y
	H <- CL_H(x, s, mu, A, b, inds)
	nrm_H <- two_norm(H)
	oldnrm_H <- nrm_H + 10*min_progress
	gamma <- .95*min(1,1/nrm_H)
	iter <- 0
	
	while(TRUE){
		
		#Step 1
		if (nrm_H < zero_tol){
			code <- 0
			cat("Solution achieved within tolerance for socp.\n")
			break
		}else if(oldnrm_H - nrm_H < min_progress){
			code <- 2
			cat("Minimum progress not achieved for socp.\n")
			break
		}else if(iter >= max_iter){
			code <- 3
			cat("Maximum number of iterations reached for socp.\n")
			break
		}else{
			rho <- CL_rho(nrm_H, gamma)
		}
		
		#Step 2
		#print(rcond(CL_grad_H(x, y, s, mu, A, c, inds), useInve = TRUE))
		#First try Newton's method
		delta_z <- try(solve(CL_grad_H(x, y, s, mu, A, c, inds), rho*zbar - H, tol=.Machine$double.xmin))
		if (class(delta_z)=="try-error"){
			code <- 1
			cat("Singularity occured for socp.")
			
			#If it didn't work then try gradient descent instead
			#delta_z <- -1*t(CL_grad_H(x, y, s, mu, A, c, inds))%*%H
		}
		
		
		#Steps 3 and 4
		lk <- -1
		done <- FALSE
		while(!done){
			lk <- lk +1
			newx <- x + (delta^lk)*delta_z[1:k]
			newy <- y + (delta^lk)*delta_z[(k+1):(k+m)]
			newmu <- mu + (delta^lk)*delta_z[(k+m+1):(k+m+n)]
			news <- c - t(A)%*%newy
			newH <- CL_H(newx, news, newmu, A, b, inds)
			newnrm_H <- two_norm(newH)
			done <- newnrm_H <= (1-sigma*(1-gamma*mu0)*(delta^lk))*nrm_H
		}
		
		x <- newx
		y <- newy
		mu <- newmu
		s <- news
		H <- newH
		oldnrm_H <- nrm_H
		nrm_H <- newnrm_H
		iter <- iter + 1
		cat("Iteration ", iter, " complete.  Magnitude = ",nrm_H,".\n",sep="")
	}
	
	
	return(list(x=x, y=y, s=s, obj=sum(c*x), code=code, mu=mu, iter=iter))
	
}

